#	File name: AJPS_Code.r

#	Author: Ryan Black 

#	Creation Date: 8/14/12

#	Last Edited: 3/30/14

##########################
#### BEGIN CODE USED #####
##########################

library(foreign)
library(cem)
library(MASS)

# Set number of values used to get L1 values. 
ND <- 1000 

# Bring in the two separate files -- one for the shortlist judges (SL) and one for the non-shortlist judges (NSL)
SL <- read.dta("~/desktop/Contender.dta", convert.factors = F, convert.underscore = F)

NSL <- read.dta("~/desktop/NonContender.dta", convert.factors = F, convert.underscore = F)

# Define outcome variables, quantiles of interest, and graphing names to be used later in making a figure
cis <- c(0.025, 0.05, 0.10, 0.50, 0.90, 0.95, 0.975)
DVs <- c("presIdeoVote", "proUS", "wrotedissent", "wroteconcur")
graph.names <- c("Presidential\nIdeological Vote",  "Pro-U.S.\nDecision", "Judge Writes\nDissent", "Judge Writes\nConcurrence")

# Create a place to store the results; the lattermost dimension is for each of the results (shortlist and non-shortlist)
gres <- array(data = NA, dim = c(length(DVs), length(cis), 7, 2))

# Do the shortlist judges first
match.vars <- c("judgeJCS", "presDist", "panelDistJCS", "circmed", "sctmed", "coarevtc", "casepub") # the set of variables to include in the matching
mods.SL <- list() # a place to store regression results
base.X <- c(0.58, 0.18, 1.17, 0.19, 0.03, 1, 1) # define baseline for pps
treat.loc <- grep("^treatFinal0$", names(SL)) # locate treatment variable
todrop <- names(SL)
keepme <- grep(paste(paste("^", match.vars, "$", sep = ""), collapse = "|"), names(SL))
todrop <- todrop[-keepme]
match.SL <- cem(treatment = "treatFinal0", drop = todrop, data = SL)
if(ND > 0){set.seed(01081982); preimb.SL <- L1.profile(data = SL, group = SL$treatFinal0, drop = todrop, plot = F, M = ND); postimb.SL <- L1.profile(data = SL, group = SL$treatFinal0, drop = todrop, plot = F, useCP = match.SL$CP, weights = match.SL$w)}
options(warn = -1)
for(D in 1:length(DVs)){
	dv.where <- grep(paste("^", DVs[D], "$", sep = ""), names(SL))
	mod.form <- as.formula(paste(DVs[D], " ~ ", "treatFinal0", " + ", paste(match.vars, collapse = " + "), sep = ""))
	mods.SL[[D]] <- glm(data = SL, mod.form, binomial(link = "logit"), weights = match.SL$w)
	set.seed(01081982)
	betas <- mvrnorm(n = 5000, mu = t(mods.SL[[D]]$coefficients), Sigma = vcov(mods.SL[[D]]))
	mfxs <- c()
	p0s <- c()
	p1s <- c()
	for(R in 1:nrow(betas)){
		mfxs[R] <- plogis(betas[R,2])-0.50
		p0s[R] <- plogis(betas[R,] %*% c(1, 0, base.X))
		p1s[R] <- plogis(betas[R,] %*% c(1, 1, base.X))
	}
	gres[D,,1,1] <- quantile(mfxs, cis)	
	gres[D,,2,1] <- quantile(p0s, cis)
	gres[D,,3,1] <- quantile(p1s, cis)	
	gres[D,,4,1] <- quantile(p1s-p0s, cis)
	bmf <- as.formula(paste(DVs[D], " ~ ", "treatFinal0",  sep = ""))
	bmf.res <- glm(data = SL, bmf, binomial(link = "logit"), weights = match.SL$w)
	set.seed(01081982)
	bmbetas <- mvrnorm(n = 5000, mu = t(bmf.res$coefficients), Sigma = vcov(bmf.res))
	p0s <- c()
	p1s <- c()
	for(R in 1:nrow(betas)){
		p1s[R] <- plogis(bmbetas[R,2]+bmbetas[R,1])
		p0s[R] <- plogis(bmbetas[R,1])
	}
	gres[D,,5,1] <- quantile(p0s, cis)	
	gres[D,,6,1] <- quantile(p1s, cis)		
	gres[D,,7,1] <- quantile(p1s-p0s, cis)			
}
options(warn = 0)

# Next do non-shortlist judges
match.vars <- c("judgeJCS", "presDist", "panelDistJCS", "circmed", "sctmed", "coarevtc") # casepub omitted b/c all Songer's are all published
mods.NSL <- list() # a place to store regression results
base.X <- c(-0.05, 0.50, 0.28, -0.02, 0.01, 1)
treat.loc <- grep("^nSLtreatFinal0$", names(NSL)) 
todrop <- names(NSL)
keepme <- grep(paste(paste("^", match.vars, "$", sep = ""), collapse = "|"), names(NSL))
todrop <- todrop[-keepme]
match.NSL <- cem(treatment = "nSLtreatFinal0", drop = todrop, data = NSL)
if(ND > 0){set.seed(01081982); preimb.NSL <- L1.profile(data = NSL, group = NSL$nSLtreatFinal0, drop = todrop, plot = F, M = ND); postimb.NSL <- L1.profile(data = NSL, group = NSL$nSLtreatFinal0, drop = todrop, plot = F, useCP = match.NSL$CP, weights = match.NSL$w)}
options(warn = -1)
for(D in 1:length(DVs)){
	dv.where <- grep(paste("^", DVs[D], "$", sep = ""), names(NSL))
	mod.form <- as.formula(paste(DVs[D], " ~ ", "nSLtreatFinal0", " + ", paste(match.vars, collapse = " + "), sep = ""))
	mods.NSL[[D]] <- glm(data = NSL, mod.form, binomial(link = "logit"), weights = match.NSL$w)
	set.seed(01081982)
	betas <- mvrnorm(n = 5000, mu = t(mods.NSL[[D]]$coefficients), Sigma = vcov(mods.NSL[[D]]))
	mfxs <- c()
	p0s <- c()
	p1s <- c()
	for(R in 1:nrow(betas)){
		mfxs[R] <- plogis(betas[R,2])-0.50
		p0s[R] <- plogis(betas[R,] %*% c(1, 0, base.X))
		p1s[R] <- plogis(betas[R,] %*% c(1, 1, base.X))
	}
	gres[D,,1,2] <- quantile(mfxs, cis)	
	gres[D,,2,2] <- quantile(p0s, cis)
	gres[D,,3,2] <- quantile(p1s, cis)	
	gres[D,,4,2] <- quantile(p1s-p0s, cis)	
	bmf <- as.formula(paste(DVs[D], " ~ ", "nSLtreatFinal0",  sep = ""))
	bmf.res <- glm(data = NSL, bmf, binomial(link = "logit"), weights = match.NSL$w)
	set.seed(01081982)
	bmbetas <- mvrnorm(n = 5000, mu = t(bmf.res$coefficients), Sigma = vcov(bmf.res))
	p0s <- c()
	p1s <- c()
	for(R in 1:nrow(betas)){
		p1s[R] <- plogis(bmbetas[R,2]+bmbetas[R,1])
		p0s[R] <- plogis(bmbetas[R,1])
	}
	gres[D,,5,2] <- quantile(p0s, cis)	
	gres[D,,6,2] <- quantile(p1s, cis)			
	gres[D,,7,2] <- quantile(p1s-p0s, cis)			
}
options(warn = 0)
dn <- list()
dn[[1]] <- DVs
dn[[2]] <- cis
dn[[3]] <- c("mfx", "p0", "p1", "pdiff", "bp0", "bp1", "bpdiff")
dn[[4]] <- c("treatFinal0", "nSLtreatFinal0")
dimnames(gres) <- dn

# Show L1 values
preimb.SL$medianL1
postimb.SL$medianL1
1-(postimb.SL$medianL1/preimb.SL$medianL1)
preimb.NSL$medianL1
postimb.NSL$medianL1
1-(postimb.NSL$medianL1/preimb.NSL$medianL1)

# Make a predicted probability plot with both SL/NSL
PTT <- c(21, 22)
YB <- 0.1
pdf(file = "~/desktop/pp.pdf", height = 7, width = 7)
par(mar = c(4,8,1,1))
plot(c(min(gres[,4,c(2,3),])-0.02, max(gres[,4,c(2,3),])+0.08), c(1-YB, (nrow(gres)*2)+YB), type = "n", xlab = "Predicted Probability of Behavior", ylab = "", axes = F, main = )
abline(h = seq(2.5, 12.5, 2))
abline(h = 1:15, lty=3)
ypts <- list()
ypts[[1]] <- rev(seq(2, nrow(gres)*2, 2))
ypts[[2]] <- rev(seq(1, nrow(gres)*2, 2))
ypts[[3]] <- rev(seq(1.5, nrow(gres)*2, 2))
for(D in 1:4){
	for(V in 1:2){
		points(gres[D,4,2,V], ypts[[V]][D], pch = PTT[V], cex = 1.69, bg = "white")
		text(gres[D,4,2,V], ypts[[V]][D]+0.2, format(round(gres[D,4,2,V], 2), nsmall = 2)) # annotate plot with point estimate values
		points(gres[D,4,3,V], ypts[[V]][D], pch = PTT[V], cex = 1.69, bg = "black")
		text(gres[D,4,3,V], ypts[[V]][D]-0.2, format(round(gres[D,4,3,V], 2), nsmall = 2)) # annotate plot with point estimate values
	}
}
axis(1, at = seq(-1,1,by=0.1))
axis(2, at = ypts[[3]], labels = graph.names, las = 2, tick = F)
box()
legend("bottomright", pch = c(19, 1, 15, 0), c("Contender, Vacancy", "Contender, No Vacancy", "Non-Contender, Vacancy", "Non-Contender, No Vacancy"), pt.cex = 1.1, bg = "white")
text(0.20, 7.1, "Difference Not Significant", font = 3)
text(0.45, 5.1, "Difference Not Significant", font = 3)
text(0.25, 3.1, "Difference Not Significant", font = 3)
text(0.45, 2.1, "Difference Not Significant", font = 3)
text(0.17, 1.1, "Difference Not Significant", font = 3)
YV <- c(4,6,8)
GL <- c(3,2,1)
for(V in 1:length(YV)){
	pc <- round(((gres[GL[V],4,3,1]/gres[GL[V],4,2,1])-1)*100,0)
	text(max(gres[,4,c(2,3),])+0.05, YV[V]+0.1, paste("+", pc, "%", sep = ""), font = 2)
}
dev.off()

########################
#### END CODE USED #####
########################
